﻿using System;
using System.Numerics;

namespace HydrogenAtom4 {
    class RealPolynomialOneHalfCoeffPrecision {
        /// <summary>
        /// 0th element: constant
        /// 1st element: 1/2th order coeff x^(1/2)
        /// 2nd element: 1st   order coeff x
        /// 3rd element: 3/2th order coeff x^(3/2)
        /// 4th element: 2nd   order coeff x^2
        ///  …
        /// </summary>
        private readonly double[] mCoefficients;

        public RealPolynomialOneHalfCoeffPrecision(double[] coeffs) {
            mCoefficients = coeffs.Clone() as double[];
        }

        public RealPolynomialOneHalfCoeffPrecision MulX() {
            // mul x means coefficients slide 2 elements

            var coeffs = new double[mCoefficients.Length+2];
            for (int i=0; i<mCoefficients.Length; ++i) {
                coeffs[2+i] = mCoefficients[i];
            }
            return new RealPolynomialOneHalfCoeffPrecision(coeffs);
        }

        public RealPolynomialOneHalfCoeffPrecision Mul(double k) {
            var coeffs = mCoefficients.Clone() as double[];
            for (int i=0; i<coeffs.Length; ++i) {
                coeffs[i] *= k;
            }
            return new RealPolynomialOneHalfCoeffPrecision(coeffs);
        }

        public double[] Coefficients {
            get { return mCoefficients; }
        }

        public RealPolynomialOneHalfCoeffPrecision Add(RealPolynomialOneHalfCoeffPrecision rhs) {
            int order = (mCoefficients.Length < rhs.mCoefficients.Length)
                ? rhs.mCoefficients.Length : mCoefficients.Length;

            var coeffs = new double[order];
            for (int i=0; i < order; ++i) {
                double v = 0.0;
                if (i < mCoefficients.Length) {
                    v += mCoefficients[i];
                }
                if (i < rhs.mCoefficients.Length) {
                    v += rhs.mCoefficients[i];
                }
                coeffs[i] = v;
            }
            return new RealPolynomialOneHalfCoeffPrecision(coeffs);
        }

        public RealPolynomialOneHalfCoeffPrecision Sub(RealPolynomialOneHalfCoeffPrecision rhs) {
            int order = (mCoefficients.Length < rhs.mCoefficients.Length)
                ? rhs.mCoefficients.Length : mCoefficients.Length;

            var coeffs = new double[order];
            for (int i=0; i < order; ++i) {
                double v = 0.0;
                if (i < mCoefficients.Length) {
                    v += mCoefficients[i];
                }
                if (i < rhs.mCoefficients.Length) {
                    v -= rhs.mCoefficients[i];
                }
                coeffs[i] = v;
            }
            return new RealPolynomialOneHalfCoeffPrecision(coeffs);
        }

        public RealPolynomialOneHalfCoeffPrecision Differential() {
            if (2 <= mCoefficients.Length && !Util.AlmostZero(mCoefficients[1])) {
                // does not support differential of 1/2th order element
                throw new NotSupportedException();
            }
            if (mCoefficients.Length <= 1) {
                // differential of constant = 0
                return new RealPolynomialOneHalfCoeffPrecision(new double[] { 0.0f });
            }
            var coeffs = new double[mCoefficients.Length - 2];
            for (int i=0; i < coeffs.Length; ++i) {
                // i=0: x       → 1             : coeffs[0] = (2/2) * mCoeffs[2]
                // i=1: x^(3/2) → (3/2)x^(1/2)  : coeffs[1] = (3/2) * mCoeffs[3]
                // i=2: x^2     → 2x            : coeffs[2] = (4/2) * mCoeffs[4]
                // i=3: x^(5/2) → (5/2)x^(3/2)  : coeffs[3] = (5/2) * mCoeffs[5]
                // i=4: x^3     → 3x^2          : coeffs[4] = (6/2) * mCoeffs[6]
                coeffs[i] = (i+2.0)/2.0 * mCoefficients[i + 2];
            }
            return new RealPolynomialOneHalfCoeffPrecision(coeffs);
        }

        public RealPolynomialOneHalfCoeffPrecision NthDifferential(int nth) {
            if (nth < 0) {
                throw new ArgumentException("nth");
            }

            var p = this;
            for (int i=0; i < nth; ++i) {
                p = p.Differential();
            }
            return p;
        }

        public Complex ValueOf(double x) {
            Complex v = mCoefficients[0];
            for (int i=1; i < mCoefficients.Length; ++i) {
                v += Util.PowComplex(x, i, 2) * mCoefficients[i];
            }
            return v;
        }
    }
}
